function [TrainTMatrix,TrainEPara,logliks] = MHMMTRAIN(Dyn,GuessTMatrix,GuessEPara,FixTMatrix,FixEPara,fhandle,varargin)

tol = 1e-6;
trtol = tol;
etol = tol;
maxiter = 500;

[numStates, checkTr] = size(GuessTMatrix);
if checkTr ~= numStates
    error('stats:hmmtrain:BadTransitions','TRANSITION matrix must be square.');
end

% number of rows of e must be same as number of states

[checkE, numParas] = size(GuessEPara);
if checkE ~= numStates
    error('stats:hmmtrain:InputSizeMismatch',...
        'EMISSIONS matrix must have the same number of rows as TRANSITIONS.');
end
if (numStates ==0 || numParas == 0)
    GuessTMatrix = [];
    GuessEPara = [];
    return
end

% initialize the counters
%TMatrix = zeros(size(GuessTMatrix));
%EPara = zeros(numStates,numParas);
converged = false;

loglik = 1; % loglik is the log likelihood of all sequences given the TR and E
logliks = zeros(1,maxiter);

oldGuessEPara = GuessEPara;
oldGuessTMatrix = GuessTMatrix;

for iteration = 1:maxiter
    oldLL = loglik;
    loglik = 0;
    
    oldGuessEPara = GuessEPara;
    oldGuessTMatrix = GuessTMatrix;

    [estimatedValues, estimatedStates,logP]  = MHMMVITERBI(Dyn,GuessTMatrix,GuessEPara,fhandle);
    loglik = loglik + logP;
    % w = warning('off');
    [TMatrix, EPara] = MHMMESTIMATE(Dyn,estimatedValues);
    % warning(w);
    % deal with any possible NaN values
    TMatrix(isnan(TMatrix)) = 0;
    totalTransitions = sum(TMatrix,2);
    % avoid divide by zero warnings
    GuessTMatrix  = TMatrix./(repmat(totalTransitions,1,numStates));
    GuessEPara = EPara;
     
    %FIX
    GuessEPara = oldGuessEPara.*FixEPara + GuessEPara.*(1-FixEPara);
    GuessTMatrix = oldGuessTMatrix.*FixTMatrix + GuessTMatrix.*(1-FixTMatrix);
    
    % Durbin et al recommend loglik as the convergence criteria  -- we also
    % use change in TR and E. Use (undocumented) option trtol and
    % etol to set the convergence tolerance for these independently.
    %
    logliks(iteration) = loglik;
    if (abs(loglik-oldLL)./(1+abs(oldLL))) < tol
        if norm(GuessTMatrix - oldGuessTMatrix,inf)./numStates < trtol
            if norm(GuessEPara - oldGuessEPara,inf)./numParas < etol
                %fprintf('Algorithm converged after %d iterations.',iteration);
                converged = true;
                break
            end
        end
    end
end

TrainTMatrix = oldGuessTMatrix.*FixTMatrix + GuessTMatrix.*(1-FixTMatrix);
TrainEPara = oldGuessEPara.*FixEPara + GuessEPara.*(1-FixEPara);;

if ~converged
    warning('stats:hmmtrain:NoConvergence',...
        'Algorithm did not converge with tolerance %f in %d iterations.',tol,maxiter);
end
logliks(logliks ==0) = [];
